Identification of Exosome-Related Genes Associated with Prognosis and Immune Infiltration Features in Head-Neck Squamous Cell Carcinoma

The highly immunosuppressive nature of head–neck squamous cell cancer (HNSCC) is not fully understood. Exosomes play crucial roles in the communication between cancer and non-cancer cells, but the clinical significance of the expression of exosome-related genes (ERGs) remains unclear in HNSCC. This study aimed to establish an HNSCC-ERGs model by using mass spectrometry (MS)-based label-free quantitative proteomics in combination with the TCGA primary HNSCC dataset. The study managed to classify the HNSCC patients into two subtypes based on the expression level of prognostic ERGs, which showed significant differences in prognosis and immune infiltration. LASSO regression algorithm was used to establish a risk prediction model based on seven risky genes (PYGL, ACTN2, TSPAN15, EXT2, PLAU, ITGA5), and the high-risk group was associated with poor survival prognosis and suppressive immune status. HPRT1 and PYGL were found to be independent prognostic factors through univariate and multivariate Cox regression analyses. Immune and ssGSEA analysis revealed that HPRT1 and PYGL were significantly related to immunosuppression, immune response, and critical signaling transduction pathways in HNSCC. Immunohistochemistry results further validated the expression level, clinical value, and immunosuppressive function of HPRT1 and PYGL in HNSCC patients. In conclusion, this study established molecular subtypes and a prediction risk model based on the ERGs. Furthermore, the findings suggested that HPRT1 and PYGL might play critical roles in reshaping the tumor microenvironment.


Introduction
The tumor microenvironment (TME) is a crucial factor in tumor progression and response to immunotherapy [1,2]. Targeting the related biological process may represent a promising approach to improving anti-tumor therapy, as suggested by Herrin [3]. Despite recent improvements in the treatment of head and neck squamous cell carcinoma (HNSCC), this disease remains highly lethal, with an overall 5-year survival rate of approximately 50% [4]. Immunotherapy has shown limited efficacy, benefiting only 15% of HNSCC patients, highlighting the need for a deeper understanding of the immunosuppressive mechanisms involved [5]. Although there have been some studies investigating the immune evasion mechanisms of cancer cells, our understanding of the maintenance of an immunosuppressive microenvironment remains incomplete and requires further investigation. R packages "ClusterProfiler" and "GSVA," employing the ssGSEA method and Spearman correlations to compute several functional pathways.

Identification of Prognostic HNSCC Exosome-Related Genes
In order to ascertain the prognostic significance of the HNSCC exosome-related genes, we utilized Cox regression analysis to assess the associations between the expression of each gene and the survival outcomes based on the TCGA-HNSCC cohort. We performed log-rank tests and summarized p-values and hazard ratios (HRs) with 95% confidence intervals (CIs). A gene was considered to have prognostic value if its p-value was less than 0.05, and we identified a set of 25 genes for further analysis.

Exploration of Gene Alteration Landscape of Prognostic ERGs
The cBio Cancer Genomics Portal (cBioPortal, https://www.cbioportal.org/, accessed on 20 September 2022) is a web portal for visual analysis of genomics in the TCGA database [21]. In this study, the cBioPortal platform was used to analyze the the effects of genetic alterations in the 25 identified prognostic ERGs on the prognosis of HNSCC patients, in terms of overall survival (OS) and disease-specific survival (DSS). The statistical significance of the results was assessed with a p-value threshold of <0.05.

Construction of Molecular Subgroups Using Consensus Clustering
The study utilized the "ConsensusClusterPlus" package to cluster the 25 prognostic ERGs genes into different subtypes in 504 HNSCC samples obtained from the TCGA database. The optimal k value was used to split the TCGA-HNSCC cohort into two clusters: cluster 1 (ERGs low expression group) and cluster 2 (ERGs high expression group). Principal component analysis (PCA) was conducted to discern the differences between the two patient groups. The survival difference between clusters 1 and 2 was evaluated using Kaplan-Meier plots.

Immune Infiltration Analysis
In order to evaluate the immune infiltration results across different subgroups, we utilized the "immuneeconv" R software package [22], which encompasses a variety of state-of-the-art algorithms such as TIMER [23], MCPCOUNTER [24], CIBERSORT [25], EPIC [26], and QUANTISEQ [27]. The immune analysis outcomes were visualized through the use of R packages "ggplot2" and "pheatmap."

Construction of Risk Model and Nomograms Based on Prognostic HNSCC ERGs
The current study aimed to investigate the potential association between the expression of HNSCC ERGs and prognosis in HNSCC patients through univariate Cox regression analysis (p < 0.05). Furthermore, a risk model was developed for HNSCC ERGs using the least absolute shrinkage and selection operator (LASSO) regression algorithm. The m6A/m1A/m5C-regulated genes were utilized to calculate the risk score (risk score = S (ExpmRNAn × βmRNAn)) [28]. The high-and low-risk groups were then compared using Kaplan-Meier curves and p-values were generated using log-rank tests. A forest plot was employed to display the p-value, HR, and 95% CI of each variable using the "forestplot" R package. Finally, a nomogram was created based on the results of multivariate Cox proportional hazards analysis of the HNSCC ERGs.

Immunohistochemistry (IHC) Staining
The IHC staining protocol was executed with meticulous adherence to our previously reported methodology [9]. Post-deparaffinization, rehydration, and antigen retrieval steps were performed, followed by quenching of endogenous peroxidase activity. Protein expression was determined using primary antibodies (rabbit anti-human HPRT1 antibody (#15059-1-AP), rabbit anti-human PYGL antibody (#15851-1-AP), and mouse anti-human CD8 antibody, all of which were acquired from Proteintech (Wuhan, China). To assess the extent of immunoreactivity, the immunoreaction score (IRS) was computed by multiplying the percentage of positive cells and the staining intensity, as previously reported [29].

Identification of HNSCC Exosome-Related Genes
Exosomes have emerged as significant regulators of intercellular communication. The present study aimed to establish a set of exosome-related genes in HNSCC. Exosomes were isolated from the cell culture media of three HNSCC cell lines, namely CAL27, SCC9, and SCC25. These genes were then overlapped with the differentially expressed genes (DEGs) of the TCGA-HNSCC cohort for further analysis ( Figure 1A). A total of 636 proteins were identified that were shared in exosomes derived from the CAL27, SCC9, and SCC25 cell lines ( Figure 1B). Consequently, the genes encoding these proteins were matched using the National Center for Biotechnology Information (NCBI) database. DEGs were identified from the TCGA database at a level of fold change > |2| and p-value < 0.05, leading to the identification of 2172 genes ( Figure 1C, Table S1). A Venn diagram was used to determine the intersection between the exosome genes from the HNSCC cell lines and DEGs from the TCGA-HNSCC cohort, which revealed 120 overlapping genes ( Figure 1D). The expression levels of these 120 HNSCC ERGs are presented in Table S2.

Functional Enrichment Analysis and Determination of Prognostic HNSCC ERGs
To ascertain the functions of the identified HNSCC ERGs, a comprehensive analysis was carried out via Gene Ontology (GO) and KEGG functional enrichment analyses. The GO analysis results showed that the 120 HNSCC ERGs were primarily implicated in cellular adhesion, migration, differentiation of epidermal/epithelial cells, response to hypoxia, extracellular exosome, and protein binding (Figure 2A-C). In parallel, the KEGG

Functional Enrichment Analysis and Determination of Prognostic HNSCC ERGs
To ascertain the functions of the identified HNSCC ERGs, a comprehensive analysis was carried out via Gene Ontology (GO) and KEGG functional enrichment analyses. The GO analysis results showed that the 120 HNSCC ERGs were primarily implicated in cellular adhesion, migration, differentiation of epidermal/epithelial cells, response to hypoxia, extracellular exosome, and protein binding (Figure 2A-C). In parallel, the KEGG pathway analysis results indicated that the 120 HNSCC ERGs were prominently associated with focal adhesion, HPV infection, the PI3K/AKT pathway, ECM-receptor interaction, and pathways in cancer ( Figure 2D). To identify HNSCC ERGs with prognostic implications, the TCGA-HNSCC cohort was subjected to univariate Cox regression analysis, which revealed 1713 genes as potential prognostic markers. Subsequently, 25 genes were identified as prognostic HNSCC ERGs associated with increased risk, with HRs > 1 ( Figure 3A,B). The prognostic value of these HNSCC ERGs was further validated by Kaplan-Meier survival curves in high-/lowexpression groups of HNSCC patients. The representative Kaplan-Meier survival curves were obtained for 9 genes ( Figure 3C-K). To identify HNSCC ERGs with prognostic implications, the TCGA-HNSCC cohort was subjected to univariate Cox regression analysis, which revealed 1713 genes as potential prognostic markers. Subsequently, 25 genes were identified as prognostic HNSCC ERGs associated with increased risk, with HRs > 1 ( Figure 3A,B). The prognostic value of these HNSCC ERGs was further validated by Kaplan-Meier survival curves in high-/lowexpression groups of HNSCC patients. The representative Kaplan-Meier survival curves were obtained for 9 genes ( Figure 3C-K).

Expression Profiles of 25 Prognostic HNSCC ERGs
Based on the fold change analysis, we observed that most of the ERGs were upregulated in HNSCC tissues relative to normal tissues, except for ACTN2, which was downregulated ( Figure 4A). To further understand the genetic alterations of HNSCC ERGs, we utilized the cBioPortal website and found that 72.28% of the 523 HNSCC cases had ERG alterations ( Figure 4B). The genetic alterations included mutation (7.84%), structural Biomolecules 2023, 13, 958 6 of 20 variant (0.38%), amplification (3.82%), deep deletion (0.96%), mRNA high (33.08%), and multiple alterations (26.2%) ( Figure 4C). Our analysis revealed that these genetic alterations significantly affected the overall survival ( Figure 4D, p < 0.001) and disease-specific survival ( Figure 4E, p < 0. 001) of HNSCC patients. Additionally, we performed Spearman's correlation analysis to analyze the correlations among the 25 prognostic ERGs, and the results are shown in Figure S2.

Expression Profiles of 25 Prognostic HNSCC ERGs
Based on the fold change analysis, we observed that most of the ERGs were upregulated in HNSCC tissues relative to normal tissues, except for ACTN2, which was downregulated ( Figure 4A). To further understand the genetic alterations of HNSCC ERGs, we utilized the cBioPortal website and found that 72.28% of the 523 HNSCC cases had ERG alterations ( Figure 4B). The genetic alterations included mutation (7.84%), structural variant (0.38%), amplification (3.82%), deep deletion (0.96%), mRNA high (33.08%), and multiple alterations (26.2%) ( Figure 4C). Our analysis revealed that these genetic alterations significantly affected the overall survival ( Figure 4D, p < 0.001) and disease-specific survival ( Figure 4E, p < 0. 001) of HNSCC patients. Additionally, we performed Spearman's correlation analysis to analyze the correlations among the 25 prognostic ERGs, and the results are shown in Figure S2.

Molecular Subtype of HNSCC Based on 25 Prognostic ERGs
To examine the association between the expression levels of HNSCC exosome-related genes and molecular subtypes of HNSCC, we conducted a consensus clustering analysis using data from the TCGA cohort comprising 504 HNSCC patients ( Figure 5A-C). Herein, we managed to identify an optimal k value of 2, and the HNSCC patients were divided into two clusters: cluster 1 with low ERG expression and cluster 2 with high ERG expression ( Figure 5D). The expression patterns of ERGs across different HNSCC subtypes are shown in Figure 5E. Notably, cluster 1 demonstrated higher overall survival ( Figure 5F), progressionfree survival ( Figure 5G) and disease-specific survival ( Figure 5H) rates than cluster 2, and the differences were statistically significant (p < 0.001). Additionally, significant differences were observed between the two clusters in terms of lymph node metastasis and grade stage (p < 0.05, Table 1). The differential expression analysis revealed that cluster 2 had 522 upregulated genes and 82 downregulated genes relative to cluster 1 ( Figure S3A,B).
KEGG pathway analysis indicated that cluster 2 was associated with activation of the TGF-β signaling pathway, regulation of actin cytoskeleton, focal adhesion, the PI3K-Akt signaling pathway, ECM receptor interaction, and suppression of tyrosine metabolism, tight junction, steroid biosynthesis, glycolysis, and the estrogen signaling pathway ( Figure S3C). The Gene Ontology analysis showed that the upregulated genes were involved in extracellular matrix organization/disassembly, endodermal cell differentiation, collagen metabolic process, and cell-substrate adhesion ( Figure S3D).

Molecular Subtype of HNSCC Based on 25 Prognostic ERGs
To examine the association between the expression levels of HNSCC exosome-related genes and molecular subtypes of HNSCC, we conducted a consensus clustering analysis using data from the TCGA cohort comprising 504 HNSCC patients ( Figure 5A-C). Herein, we managed to identify an optimal k value of 2, and the HNSCC patients were divided into two clusters: cluster 1 with low ERG expression and cluster 2 with high ERG expression ( Figure 5D). The expression patterns of ERGs across different HNSCC subtypes are ure S3A,B). KEGG pathway analysis indicated that cluster 2 was associated with activat of the TGF-β signaling pathway, regulation of actin cytoskeleton, focal adhesion, the PI Akt signaling pathway, ECM receptor interaction, and suppression of tyrosine meta lism, tight junction, steroid biosynthesis, glycolysis, and the estrogen signaling pathw ( Figure S3C). The Gene Ontology analysis showed that the upregulated genes were volved in extracellular matrix organization/disassembly, endodermal cell differentiati collagen metabolic process, and cell-substrate adhesion ( Figure S3D).

HNSCC ERG-High Subtype Was Associated with Immunosuppressive TME
Exosomes play a crucial role in regulating intercellular communication. Here, we investigated the characteristics of immune infiltration in the tumor microenvironment (TME) of HNSCC. Our analysis revealed a significant difference between the ERG-high cluster (cluster 2) and ERG-low cluster (cluster 1) in terms of low immune score ( Figure 6A) and microenvironment score ( Figure 6C) but not stroma score ( Figure 6B). We utilized the TIDE algorithm to predict the tumor immune response rate and found that cluster 2 had a significantly higher TIDE score than cluster 1 (p = 1.4 × 10 −15 ), suggesting that patients in the ERG-low group might benefit more from immune checkpoint inhibitor therapy (ICBs) ( Figure 6D). Using the CIBERSORT algorithm, we analyzed immune heterogeneity between the two molecular subgroups ( Figure 6E). The results showed that cluster 2 was negatively associated with B cell plasma, CD8+ T cells, activated CD4+ memory T cells, and activated NK cells (p < 0.001) but positively associated with resting CD4+ memory T cells and resting NK cells (p < 0.001). Further validation using the MCPCOUNTER algorithm showed that cluster 2 was associated with stronger immunosuppressive status, with a positive relationship with T cells, CD8+ T cells, and B cells (p < 0.01) ( Figure 6F). Similar results were obtained using the EPIC algorithm ( Figure S4A), QUANTISEQ algorithm ( Figure S4B), and TIMER algorithm ( Figure S4C). We also analyzed the expression levels of immune checkpoint markers in both HNSCC subtypes and found that the expression of CD274, HAVCR2, PDCD1LG2, and SIGLEC15 was elevated in cluster 2 ( Figure 6G). Overall, our findings suggested that HNSC patients in the ERG-high subtype might have a greater likelihood of developing an immunosuppressive microenvironment, which could ultimately contribute to poor prognosis. results were obtained using the EPIC algorithm ( Figure S4A), QUANTISEQ algorithm ( Figure S4B), and TIMER algorithm ( Figure S4C). We also analyzed the expression levels of immune checkpoint markers in both HNSCC subtypes and found that the expression of CD274, HAVCR2, PDCD1LG2, and SIGLEC15 was elevated in cluster 2 ( Figure 6G). Overall, our findings suggested that HNSC patients in the ERG-high subtype might have a greater likelihood of developing an immunosuppressive microenvironment, which could ultimately contribute to poor prognosis.

Construction of 25 Prognostic HNSCC ERGs Risk Model
To further investigate the association between ERGs and prognosis in HNSCC, we employed the least absolute shrinkage and selection operator (LASSO) Cox regression algorithm to determine the minimum λ value ( Figure 7A,B). Subsequently, we constructed a 7-gene signature based on the optimal λ value. The correlation between the expression levels of these 7 genes and the risk score is displayed in Figure 7C. The risk score was calculated as follows: Riskscore = (0.0944) × PYGL + (0.0237) × ACTN2 + (0.0133) × TSPAN15 + (0.0161) × EXT2 + (0.0368) × PLAU + (0.0092) × ITGA5 + (0.3567) × HPRT1. Based on the risk score, HNSCC patients were stratified into two groups, and Kaplan-Meier survival analysis indicated that patients in the high-risk group had a worse prognosis than those in the low-risk group ( Figure 7D). a 7-gene signature based on the optimal λ value. The correlation between the expression levels of these 7 genes and the risk score is displayed in Figure 7C. The risk score was calculated as follows: Riskscore = (0.0944) × PYGL + (0.0237) × ACTN2 + (0.0133) × TSPAN15 + (0.0161) × EXT2 + (0.0368) × PLAU + (0.0092) × ITGA5 + (0.3567) × HPRT1. Based on the risk score, HNSCC patients were stratified into two groups, and Kaplan-Meier survival analysis indicated that patients in the high-risk group had a worse prognosis than those in the low-risk group ( Figure 7D).

Correlation between Immune Score and Risky ERGs in HNSCC
To assess the potential relationship between the 7-gene risk model and immune score in HNSCC, we utilized the MCPCOUNTER ( Figure 8) and TIMER algorithms ( Figure S5) for analysis. Our results showed that the risk score was significantly and negatively associated with the infiltration of T cells ( Figure 8A), CD8+ T cells ( Figure 8B), NK cells ( Figure  8D), B cells ( Figure 8E), and myeloid dendritic cells ( Figure 8H) in the tumor microenvironment, while no significant differences were found for the level of cytotoxic lymphocytes ( Figure 8C), monocytes ( Figure 8F), and neutrophils ( Figure 8I). Interestingly, the risk score was positively correlated with endothelial cells, suggesting that HNSCC-derived exosomes might contribute to angiogenesis ( Figure 8G). We further evaluated the

Correlation between Immune Score and Risky ERGs in HNSCC
To assess the potential relationship between the 7-gene risk model and immune score in HNSCC, we utilized the MCPCOUNTER ( Figure 8) and TIMER algorithms ( Figure S5) for analysis. Our results showed that the risk score was significantly and negatively associated with the infiltration of T cells ( Figure 8A), CD8+ T cells ( Figure 8B), NK cells ( Figure 8D), B cells ( Figure 8E), and myeloid dendritic cells ( Figure 8H) in the tumor microenvironment, while no significant differences were found for the level of cytotoxic lymphocytes ( Figure 8C), monocytes ( Figure 8F), and neutrophils ( Figure 8I). Interestingly, the risk score was positively correlated with endothelial cells, suggesting that HNSCC-derived exosomes might contribute to angiogenesis ( Figure 8G). We further evaluated the correlation between the expression level of the 7 risky ERGs and immune cell infiltration. We used 4 algorithms (MCPCOUNTER, TIMER, QUANTISEQ, and EPIC) for analysis in HNSCC and found that the expression levels of HPRT1 and PYGL were related to decreased CD8+ T cell infiltration level across all four algorithms ( Figure 9A-D). We also examined the correlation between the expression profile of the 7 risky genes and CD8+ T cell markers (CD3D, CD3E, CD8A, CD8B, IFNG, GZMB) and found that EXT2, HPRT1, ITGA5, PYGL, and PLAU were negatively associated with CD8+ T cell abundance. The correlation between the 7 risky genes and other infiltration cell markers was also explored ( Figure S6). These results suggested that the 7-gene risk model might be associated with immune cell infiltration in HNSCC, and the identified risky genes might play a potential role in modulating the immune response in HNSCC. examined the correlation between the expression profile of the 7 risky genes and CD8+ T cell markers (CD3D, CD3E, CD8A, CD8B, IFNG, GZMB) and found that EXT2, HPRT1, ITGA5, PYGL, and PLAU were negatively associated with CD8+ T cell abundance. The correlation between the 7 risky genes and other infiltration cell markers was also explored ( Figure S6). These results suggested that the 7-gene risk model might be associated with immune cell infiltration in HNSCC, and the identified risky genes might play a potential role in modulating the immune response in HNSCC.

A prognostic Nomogram Based on the 7 Risky ERGs
In our study, we developed a prognostic nomogram that utilized the expression levels of the 7 risky genes to predict the survival probability in HNSCC patients. Univariate ( Figure 10A) and multivariate ( Figure 10B) Cox analyses were conducted to assess the relationship between individual genes and patient prognosis. The results showed that HPRT1 and PYGL might be independent prognostic factors in HNSCC. Based on these findings, a nomogram was constructed that incorporated the expression levels of HPRT1 and PYGL ( Figure 10C). The predictive nomogram was able to accurately predict the 1-, 3-, and 5-year overall survival rates in the entire cohort ( Figure 10D). Together with the immune analysis presented in Figure 9, our results suggested that HPRT1 and PYGL might be potential key ERGs in HNSCC.

A prognostic Nomogram Based on the 7 Risky ERGs
In our study, we developed a prognostic nomogram that utilized the expression levels of the 7 risky genes to predict the survival probability in HNSCC patients. Univariate ( Figure 10A) and multivariate ( Figure 10B) Cox analyses were conducted to assess the relationship between individual genes and patient prognosis. The results showed that HPRT1 and PYGL might be independent prognostic factors in HNSCC. Based on these findings, a nomogram was constructed that incorporated the expression levels of HPRT1 and PYGL ( Figure 10C). The predictive nomogram was able to accurately predict the 1-, 3-, and 5-year overall survival rates in the entire cohort ( Figure 10D). Together with the immune analysis presented in Figure 9, our results suggested that HPRT1 and PYGL might be potential key ERGs in HNSCC.

Genetic Function Analysis of HPRT1 and PYGL
The present study explored the GEPIA database to investigate the expression levels and functions of HPRT1 and PYGL in multiple cancer types. HPRT1 expression was found to be upregulated in several cancer types, including BRCA, CESC, COAD, DLBC, ESCA, HNSC, LUAD, LUSC, PAAD, READ, SKCM, STAD, THYM, and UCEC, while it was downregulated in GBM, LAML, LGG, and TGCT ( Figure S7A). Similarly, PYGL expression was upregulated in GBM, HNSC, KIRC, KIRP, LAML, LGG, PAAD, SKCM, and TGCT, but it was downregulated in ACC, DLBC, KICH, and THYM ( Figure S7B). Furthermore, gene set enrichment analysis (GSEA) revealed that HPRT1 was positively associated with MYC targets, G2M checkpoint, tumor proliferation signature, DNA repair, cellular response to hypoxia, and ferroptosis, while it was negatively associated with ECM degradation, the P53 pathway, inflammatory response, and the IL10 anti-inflammatory signaling pathway ( Figure 11B). In contrast, PYGL was positively correlated with the P53 pathway, cellular response to hypoxia, TGFβ, apoptosis, EMT marker, angiogenesis, collagen formation, and MYC targets, but it was negatively correlated with DNA replication and the tumor inflammation signature ( Figure 11B). KEGG and GO functional enrichment analyses were also conducted, and the HPRT1 high expression group was associated with activation of the Wnt signaling pathway, drug metabolism, purine metabolism, ether lipid metabolism, skin development, epidermis development, and suppression of the IL17 signaling pathway compared to the HPRT1 low expression group ( Figure S8). Similarly, the PYGL high expression group was highly associated with PI3K signaling, the Hippo pathway, human papillomavirus infection, extracellular matrix organization/assembly, and epidermis development. Overall, these findings suggested that HPRT1 and PYGL might play important roles in cancer development and prognosis ( Figure S9).

Genetic Function Analysis of HPRT1 and PYGL
The present study explored the GEPIA database to investigate the expression levels and functions of HPRT1 and PYGL in multiple cancer types. HPRT1 expression was found to be upregulated in several cancer types, including BRCA, CESC, COAD, DLBC, ESCA, HNSC, LUAD, LUSC, PAAD, READ, SKCM, STAD, THYM, and UCEC, while it was downregulated in GBM, LAML, LGG, and TGCT ( Figure S7A). Similarly, PYGL expression was upregulated in GBM, HNSC, KIRC, KIRP, LAML, LGG, PAAD, SKCM, and TGCT, but it was downregulated in ACC, DLBC, KICH, and THYM ( Figure S7B). Furthermore, gene set enrichment analysis (GSEA) revealed that HPRT1 was positively associated with MYC targets, G2M checkpoint, tumor proliferation signature, DNA repair, cellular response to hypoxia, and ferroptosis, while it was negatively associated with ECM degradation, the P53 pathway, inflammatory response, and the IL10 anti-inflammatory signaling pathway ( Figure 11B). In contrast, PYGL was positively correlated with the P53 pathway, cellular response to hypoxia, TGFβ, apoptosis, EMT marker, angiogenesis, collagen formation, and MYC targets, but it was negatively correlated with DNA replication and the tumor inflammation signature ( Figure 11B). KEGG and GO functional enrichment analyses were also conducted, and the HPRT1 high expression group was associated with activation of the Wnt signaling pathway, drug metabolism, purine metabolism, ether lipid metabolism, skin development, epidermis development, and suppression of the IL17 sig- pathway, human papillomavirus infection, extracellular matrix organization/assembly, and epidermis development. Overall, these findings suggested that HPRT1 and PYGL might play important roles in cancer development and prognosis ( Figure S9).

Overexpression of HPRT1 and PYGL Were Correlated with Tumor Progression and CD8+ T Cell Infiltration in HNSCC Patients
In this study, the protein expression levels of HPRT1 and PYGL were evaluated using immunohistochemistry in a validated primary HNSCC cohort (n = 48). Representative images of low, moderate, or high levels based on IRS score are shown in Figure 12A. Significantly higher expression levels of HPRT1 and PYGL were observed in HNSCC tissue as compared to the paired normal control (p < 0.01) ( Figure 12B). No statistical difference was observed based on age ( Figure S10). Correlation analysis was performed to investigate the relationship between expression levels of PYGL and HPRT1, tumor size, and lymph node metastasis. The expression level of PYGL was found to be positively correlated with tumor size ( Figure 12C, p < 0.05) and lymph node metastasis ( Figure 12D, p < 0.05). However, only lymph node status was correlated with HPRT1 expression ( Figure 12C,D). In addition, the study validated the relationship between HPRT1, PYGL, and CD8+ T cell infiltration level (Figure 8) in the primary HNSCC cohort. A negative correlation was observed between HPRT1 and PYGL expression levels and CD8+ T cell infiltration level ( Figure  12E). These findings suggested that HPRT1 and PYGL might be associated with tumor progression and might induce an immunosuppressive microenvironment.

Overexpression of HPRT1 and PYGL Were Correlated with Tumor Progression and CD8+ T Cell Infiltration in HNSCC Patients
In this study, the protein expression levels of HPRT1 and PYGL were evaluated using immunohistochemistry in a validated primary HNSCC cohort (n = 48). Representative images of low, moderate, or high levels based on IRS score are shown in Figure 12A. Significantly higher expression levels of HPRT1 and PYGL were observed in HNSCC tissue as compared to the paired normal control (p < 0.01) ( Figure 12B). No statistical difference was observed based on age ( Figure S10). Correlation analysis was performed to investigate the relationship between expression levels of PYGL and HPRT1, tumor size, and lymph node metastasis. The expression level of PYGL was found to be positively correlated with tumor size ( Figure 12C, p < 0.05) and lymph node metastasis ( Figure 12D, p < 0.05). However, only lymph node status was correlated with HPRT1 expression ( Figure 12C,D). In addition, the study validated the relationship between HPRT1, PYGL, and CD8+ T cell infiltration level (Figure 8) in the primary HNSCC cohort. A negative correlation was observed between HPRT1 and PYGL expression levels and CD8+ T cell infiltration level ( Figure 12E). These findings suggested that HPRT1 and PYGL might be associated with tumor progression and might induce an immunosuppressive microenvironment.

Discussion
Numerous studies have demonstrated the importance of exosomes in investigating the mechanisms underlying the development, advancement, metastasis, and immune evasion of HNSCC [30][31][32]. The complexity of the TME plays a critical role in tumor progression. In addition to cancer cells, the TME comprises a variety of non-tumor components, such as immune cells, endothelial cells, fibroblasts, and cellular metabolites. Exosomes are a specific type of EV with a wide range of biological functions in mediating cellular interactions between cancer and non-tumor cells through the proteins, nucleic acids, and metabolites they transport [33]. However, the precise roles of exosome-mediated gene regulation in the progression of HNSCC remain unclear. In this study, we aimed to investigate the biological functions and prognostic values of the HNSCC ERGs.
In this study, a HNSCC ERG set was constructed using mass spectrometry-based label-free quantitative proteomics and the TCGA-HNSCC dataset. Through differential analysis and functional enrichment analysis, 120 genes were identified as HNSCC ERGs. Among them, 25 genes were found to have a significant value in overall survival. All the included HNSCC patients were classified into two subtypes, and we observed that the ERGhigh expression cluster was associated with poor prognosis and an immunosuppressive microenvironment. A risk model and nomogram were developed based on seven risky genes. The study also revealed that HPRT1 and PYGL were the core ERGs owing to their prognostic value and potential immunosuppressive function, which was further confirmed in the validated HNSCC cohort.
Exosomes play a crucial role in tumorigenesis as essential immunomodulators, particularly in remodeling the TME to induce tumor metastasis and an immunosuppressive status [34]. For example, exosomal PDL1 derived from tumor cells leads to immunosuppression and correlates with anti-PDL1 immunotherapy [11]. Exosomal CD73 derived from HNSCC tumor cells mediates phenotypic and functional changes in TAMs to induce immune tolerance [30]. However, systematic research on exosome-related genes in cancer based on big data is still lacking. In this study, we conducted a comprehensive analysis of ERGs in the HNSCC cohort. Identifying cancer subgroups based on gene expression is meaningful in clinical practice [35]. Based on the prognostic ERGs, we reclassified HNSCC patients into two subtypes, which showed significant clinical implications. A prognostic model was established based on PYGL, ACTN2, TSPAN15, EXT2, PLAU, ITGA5, and HPRT1. Univariate and multivariate Cox regression analyses further revealed that HPRT1 and PYGL might be independent factors for constructing a nomogram model. HPRT1 was reported to play a critical role in cancer progression and chemoresistance, including HNSCC. Overexpression of PYGL also predicted poor prognosis in solid tumors. Both genes were detected in HNSCC-derived exosomes and validated in primary HNSCC tissue. In the immune infiltration analysis, both genes were significantly correlated with immune inhibition. We focused on T cells, which are the main anti-tumor immune cells. Accordingly, overexpression of HPRT1 and PYGL was negatively related to T cell (especially in CD8+ T cell) infiltration level and CD8+ T cell effectors (IFNG, GZMB). Herein, the two genes were identified as key ERGs of HNSCC. The immune surveillance system functions as a tumor suppressor in the early stage of tumor development. T cell metabolism altered by oncometabolite is a novel capability of cancer cells for immune evasion. Thus, we hypothesize that the two metabolic enzymes might reprogram T cell metabolism to sustain immune inhibition status. However, more research is needed to explore the potential mechanism.
HPRT1 is recognized as a significant contributor to cancer progression and chemoresistance, particularly in HNSCC [36]. PYGL overexpression is a predictor of poor prognosis in solid tumors [37]. HPRT1 and PYGL were detected in HNSCC exosomes and verified in primary HNSCC tissue. Upon conducting immune infiltration analysis, we found that both genes exhibited a significant correlation with immune inhibition, particularly with T cells, which are the main anti-tumor immune cells. Notably, overexpression of HPRT1 and PYGL was negatively correlated with the infiltration level of T cells, particularly CD8+ T cells, and CD8+ T cell effector molecules, including IFNG and GZMB. Therefore, we identified HPRT1 and PYGL as crucial ERGs in HNSCC. The immune surveillance system plays a crucial role in the early stages of tumor development as a tumor suppressor [38]. A novel cancer cell capability for immune evasion is altering T cell metabolism through oncometabolites [39]. Hence, we postulate that the two metabolic enzymes might reprogram T cell metabolism to sustain immune inhibition status. However, further research is required to investigate the potential mechanisms involved.
In summary, we developed a set of 120 HNSCC ERGs and identified a correlation between HNSCC reclassification, poor prognosis, and an immunosuppressive microenvironment. Additionally, we established a risk model and nomogram based on seven risky genes (PYGL, ACTN2, TSPAN15, EXT2, PLAU, ITGA5, HPRT1). Through immune infiltration analysis and determination of prognostic value, HPRT1 and PYGL were identified as key regulators carried by exosomes, based on their prognostic value and correlation with immune inhibition. However, further research is necessary to determine the mechanisms by which HPRT1 and PYGL induce immune evasion.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/biom13060958/s1, Figure S1: The size and quantity of exosomes were measured using NTA; Figure S2: The correlation among 25 prognostic ERGs was assessed via Spearman's correlation analysis; Figure Figure S10: The correlation between HPRT1, PYGL expression and age; Table S1: The DEGs in TCGA-HNSCC cohort; Table S2: The expression levels of these 120 HNSCC ERGs.